home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / gersh_r < prev    next >
Encoding:
Text File  |  1995-01-16  |  1.5 KB  |  66 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Gershgorin disks.
  4.  
  5. // Syntax:      gersh ( A )
  6.  
  7. // Description:
  8.  
  9. //      gersh(A) draws the Gershgorin disks for the matrix A.  The
  10. //      eigenvalues are plotted as crosses `x'.
  11.  
  12. //      Alternative usage: gersh(A, 1) suppresses the plot and returns
  13. //      the data in G, with A's eigenvalues in E.
  14. //
  15. //      Try: gersh(lesp(n)) and gersh(smoke(n,1)).
  16.  
  17. //    This file is a translation of gersh.m from version 2.0 of
  18. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  19. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  20.  
  21. // Dependencies
  22.    require seqa cpltaxes
  23.  
  24. //-------------------------------------------------------------------//
  25.  
  26. gersh = function (A, noplot)
  27. {
  28.   local (A, noplot)
  29.   global (pi)
  30.  
  31.   if (diff (size (A))) { error ("Matrix must be square."); }
  32.  
  33.   n = max (size (A));
  34.   m = 40;
  35.   G = zeros (m,n);
  36.  
  37.   d = diag (A);
  38.   r = sum ( abs ( A.'-diag (d) ) )';
  39.   e = eig (A).val.';
  40.  
  41.   radvec = exp (1i * seqa (0, 2*pi,m)');
  42.  
  43.   for (j in 1:n)
  44.   {
  45.     G[;j] = d[j]*ones (size (radvec)) + r[j]*radvec;
  46.   }
  47.  
  48.   if (!exist (noplot))
  49.   {
  50.     ax = cpltaxes (G[:]);
  51.     plimits (ax[1], ax[2], ax[3], ax[4]);
  52.     plaxis (); plimits (); plaspect (1);
  53.     plstyle ("line");
  54.     for (j in 1:n)
  55.     {
  56.       plimits (ax[1], ax[2], ax[3], ax[4]);
  57.       plot ([ real (G[;j]), imag(G[;j]) ]);      // Plot the disks.
  58.       subplot (0);
  59.     }
  60.     plstyle ("point");
  61.     plot ([ real (e), imag (e) ]);     // Plot the eigenvalues too.
  62.   }
  63.  
  64.   return << G=G; e=e >>;
  65. };
  66.